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ABSTRACT 


V 

The  performance  of  a  single-stage  surface -to -orbit  shuttle — 
whether  chemical -propellant  or  nuclear— can  be  considerably  improved  by 
*mixed-modet  propulsion.  A  mixed-mode  shuttle  would  be  fitted  with 
engines  designed  to  use  two  different  propellant  combinations— a  high- 
thrust  (mode  1)  propellant,  such  as  ammonia,  and  a  high-specific-impulse 
(mode  2)  propellant,  such  as  hydrogen. 

The  first  step  in  the  evaluation  of  the  mixed-mode  nuclear  shuttle 
is  a  preliminary  trajectory  optimization  study.  This  study,  using  a 
simple  mission  and  a  simple  shuttle  model,  would  be  the  basis.for  more 
complex  trajectory  optimization  studies.  The  problem  considered  in  this 
thesis  is  that  of  minimizing  the  propellant  expenditure  of  a  mixed-mode 
nuclear  shuttle  for  a  given  orbit.  The  starting  point  for  any  optimiza¬ 
tion  problem  is  a  mathematical  model  of  the  system,  in  state  variable 
form. 

Application  of  the  methods  of  optimal  control  theory  results  in  a 
two -point  boundary-value  problem  and  an  associated  algebraic  problem. 
Three  numerical  methods  are  briefly  described  which  could  be  used  to 
solve  these  problems.^  However  a  more  flexible  method,  that  of  finite 
differences,  is  proposed  in  this  study  to  handle  the  more  restricted 
problem.  Replacement  of  the  derivatives  with  finite-difference  approxi¬ 
mations  results  in  a  set  of  nonlinear  algebraic  equations.  These 
equations  are  solved  by  the  Newton-Raphson  technique;  central -difference 
approximations  are  used  for  the  Jacobian  matrix,  and  the  linearized 
algebraic  equations  are  solved  by  the  Gauss  elimination  method.  A 
procoduro  is  prosentod  for  obtaining  good  starting  values  for  the 
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Newton-Raphson  iteration.  This  is  done  by  neglecting  the  effects  of 
atmosphere;  as  a  result,  the  optimality  conditions  are  greatly 
simplified.  ,A$  it  turned  out,  this  procedure  proved  sufficient  for  an 
evaluation  of  nixed -mode  propulsion.  Consequently,  the  finite- 
difference  method  was  not  used  to  obtain  numerical  results,  although  the 
correctness  of  the  program  was  verified. 

The  analysis  of  the  problem  is  given  in  detail,  along  with  numerical 
results  for  various  combinations  of  the  input  variables,  and  the 
finite-difference  computer  program,  properly  documented. 

In  conclusion,  suggestions  for  extensions  of  this  work  are  presented 
and  methods  of  approach  are  outlined. 
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CHAPTER  I 


INTRODUCTION 

Since  the  end  of  the  Apollo  program,  the  trend  in  space  transpor¬ 
tation  has  been  towards  low  cost.  The  intention  was  that  reusability 
of  space  vehicles  and  equipment,  conaonality,  and  integration  of  manned 
and  unmanned  flights  would  accomplish  this  task.  The  space  shuttle  was 
born  out  of  this  trend. 

The  present  shuttle  was  supposed  to  drastically  reduce  the  cost  of 
surface -to -orbit  transportation.  However,  this  shuttle  will  be  cheaper 
to  develop  than  to  operate.  The  estimated  payload  costs  were  350  dollars 
per  kilogram  per  low  altitude  low  inclination  orbits  and  the  actual 
figure  will  undoubtedly  be  substantially  higher. 

The  present  shuttle  has  operational  deficiencies  as  well.  One  of 
its  main  objectives  is  to  carry  large  payloads  into  earth  orbit  with 
greater  reliability  and  economic  payback.  But  the  payload  capability 
itself  is  quite  limited— only  29,500  kilograms  maximum  in  low  altitude 
low  inclination  orbit  and  this  decreases  rapidly  as  orbital  altitude  and 
inclination  increase.  In  addition,  there  is  very  little  orbital  maneuver 
capability.  Finally,  one  has  the  nuisance  of  sea-recovering  the  boosters. 

Therefore  it  is  logical  to  look  to  the  possibility  of  a  single- 
stage  shuttle  because  operational  and  developmental  costs  could  be 
greatly  reduced.  But  the  present  propulsion  systems  do  not  permit 
single-stage  winged  shuttles  with  reasonably  large  payload  fractions  or 
perhaps  with  any  payload  at  all. 

The  concept  of  mixed-mode  propulsion  may  make  a  reusable  one-stage- 
to-orbit  shuttle  feasible.  "Mixed -mode"  propulsion  for  surface-to-orbit 
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shuttles  was  first  proposed  by  Robert  Salkeld  in  an  article  published  in 
1971.*  There  is  a  tradeoff  between  thrust  and  specific  impulse ,  but 
this  tradeoff  shifts  during  ascent.  High  thrust  is  at  a  premium  at 
liftoff,  but  as  speed  and  altitude  are  built  up  the  tradeoff  shifts  in 
the  direction  of  high  specific  impulse.  A  mixed-mode  shuttle  would  be 
fitted  with  engines  designed  to  bum  two  different  propellant  combina¬ 
tions — a  high-thrust  (mode  1)  combination,  such  as  kerosene  and  oxygen, 
and  a  high-specific-impulse  (mode  2)  combination,  such  as  hydrogen  and 
oxygen— simultaneously.  At  liftoff,  the  engines  would  be  operating  in 
mode  1;  they  would  transition  to  mode  2  during  ascent.  Detailed  studies 
have  shown  that  mixed-mode  propulsion  makes  single-stage  winged  shuttles 
of  modest  size  practical. 

Based  on  initial  economic  analysis,  the  initial  cost  of  mixed-mode 
shuttles  will  be  far  less  than  for  the  single-mode  alternatives.  Robert 
Salkeld  believes  the  mixed-mode  concept  can  cut  the  overall  cost  of  a 
1000-flight  shuttle  program  at  least  40%  while  producing  a  much  more 
effective  system. 

Operationally,  the  mixed-mode  shuttle  also  seems  to  be  a  good  bet. 
The  mixed-mode  propulsion  can  increase  performance  in  a  vehicle  of  a 
given  size  or  reduce  the  size  of  a  vehicle  of  a  given  performance. 
Propellant  volume  is  reduced,  allowing  for  a  larger  payload  capability, 
thus  satisfying  one  of  the  main  objectives  for  the  space  shuttle.  Thus 
the  potential  performance  gains  from  mixed-mode  propulsion  warrant  its 
study. 

Donald  Kingsbury  has  pointed  out  that  the  mixed-mode  concept 
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can  also  be  applied  to  nuclear  rockets.  The  most  likely  mode  1  pro¬ 
pellant  is  ammonia,  which  has  a  fairly  high  density  but  readily 


dissociates  into  a  mixture  of  nitrogen  and  hydrogen  having  a  fairly  low 
molecular  weight.  The  mode  2  propellant  would  be  hydrogen. 

The  first  step  in  the  evaluation  of  the  performance  of  a  mixed-mode 
nuclear  shuttle  is  an  ascent -trajectory  optimization  study.  The  purpose 
of  this  work  is  to  program  the  thrust  vector  for  a  mixed-mode  single- 
stage  nuclear  surface-to-orbit  shuttle  so  as  to  minimize  the  propellant 
expenditure  required  to  reach  a  given  orbit.  (Minimizing  the  propellant 
expenditure  required  to  reach  a  given  orbit  is  equivalent  to  maximizing 
the  final  mass  of  the  shuttle.)  The  following  assumptions  are  made: 

1.  The  shuttle  is  launched  from  a  site  on  the  equator  into 
an  equatorial  orbit.  As  is  appropriate  for  a  preliminary 
study,  this  assumption  is  restrictive,  but  not  unrealistic.  It 
considerably  simplifies  the  computation  and  the  presentation  of 
results,  without  sacrificing  any  of  the  essential  features  of 
the  problem.  The  shuttle's  performance  would  be  slightly  reduced 
for  non -equatorial  launch  sites  and  non -equatorial  orbits. 

2.  The  nozzle  flow  may  be  treated  as  a  one-dimensional  ideal -gas 
flow.  This  assumption  is  acceptable  for  a  preliminary  study, 
and  it  greatly  simplifies  the  computation. 

3.  The  nozzle  flow  is  expanded  to  ambient  pressure.  This  assumption 
implies  that  the  shuttle  is  fitted  with  an  aerospike  nozzle.  If 
this  is  not  the  case,  the  assumption  is  slightly  optimistic. 

4.  The  effective  ptnp  power  and  the  reactor-exit  stagnation  temper¬ 
ature  are  independent  of  the  reactor-exit  stagnation  pressure. 
("Effective  pump  power"  is  the  power  required  to  pump  the 
propellants,  assumed  to  be  incompressible,  from  zero  pressure 

to  the  reactor -exit  stagnation  pressure.)  The  accuracy  of  this 


assumption  depends  upon  the  quality  of  the  pump  and  heat* 
exchanger  designs. 

5.  The  effect  of  Reynolds  number  upon  the  drag  and  lift  coeffi¬ 
cients  is  negligible. 

6.  Atmospheric  temperatures  and  pressures  are  approximately  those 
of  the  Standard  Tropical  Atmosphere. 

The  following  axe  given,  in  addition  to  the  configuration  and  the 
choice  of  propellants: 

1.  Final  altitude  and  velocity. 

2.  Reactor-exit  stagnation  temperature. 

3.  Effective  pump  power/initial  mass. 

4.  Nozzle  throat  area/ initial  mass. 

5.  Base  area/ initial  mass. 

(Note  that,  for  a  given  family  of  designs,  the  payload/initial  mass  ratio 
may  depend  upon  all  these  things,  as  well  as  upon  the  propellant  capacity 
per  unit  initial  mass.) 

The  following  is  an  outline  of  the  remainder  of  this  work: 

CHAPTER  II.  THE  MATHEMATICAL  FORMULATION  OF  TOE  PROBLEM 

The  differential  equations  (equations  of  motion)  which  describe  the 

system  are  derived  and  explained.  Boundary  conditions  are  given. 

The  performance  index  is  defined. 

CHAPTER  III.  THE  OPTIMAL  CONTROL  PROBLEM 

Optimal  control  theory  is  applied  to  obtain  the  equations  which  must 
be  solved  to  compute  the  optimal  trajectory. 

CHAPTER  IV.  THE  NUMERICAL  SOLUTION 

A  finite-difference  method  for  the  solution  of  the  optimal  control 
problem  is  described.  Three  alternative  numerical  methods  are  also 


CHAPTER  V.  APPROXIMATE  COMPUTATION  OF  TOE  OPTIMAL  TRAJECTORY 
TOe  algorithm  for  obtaining  good  starting  values  to  be  used  to  begin 
the  Newton -Raphson  iteration. 

CHAPTER  VI .  NUMERICAL  RESULTS  AND  CONCLUSIONS 

The  results  are  presented  from  the  approximate  trajectory  method. 

CHAPTER  VII.  RECOMMENDATIONS 

Possible  extensions  of  this  work  and  methods  of  approach  are  sug- 
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CHAPTER  II 

THE  MATHEMATICAL  FORMULATION  OF  THE  PROBLEM 

The  starting  point  for  an  optimal  control  investigation  is  a  mathe¬ 
matical  model  in  state  variable  form.  A  state  vector  for  a  system  is  a 
minimum  set  of  data  required  to  predict  the  future  behavior  of  the  system. 
The  state  variables  satisfy  a  set  of  first-order  ordinary  differential 
equations  called  the  state  differential  equations.  The  state  differen¬ 
tial  equations  contain  certain  variables  in  addition  to  the  state  vari¬ 
ables.  These  are  referred  to  as  the  control  variables.  The  control 
variables  are  usually  taken  to  be  quantities  which  can  actually  be  mani¬ 
pulated. 

The  choice  of  what  variables  are  to  be  treated  as  state  variables 
and  what  variables  are  to  be  treated  as  control  variables  is  purely  one 
of  convenience.  However,  the  number  of  state  variables  and  the  number 
of  control  variables  is  determined  by  the  mathematical  model  of  the 
system. 

The  control  optimization  problem  is  that  of  finding  a  control  vector 
that  "optimizes"  the  performance  of  the  system. 

For  this  case,  let  the  state  variables  be  r,  r,  and  ♦  ,  where  r  is 

the  distance  from  the  center  of  the  planet  and  +  is  longitude;  let  the 

control  variables  be  0  and  p  .  where  0  is  the  angle  from  the  horizontal 

r 

to  the  thrust  axis,  positive  upward,  and  pr  is  the  stagnation  pressure  at 
the  reactor  exit;  and  let  the  independent  variable  be  m,  where  m  is  the 
ratio  of  the  shuttle's  instantaneous  mass  to  the  initial  mass. 

The  shuttle's  radial  and  transverse  accelerations  in  terms  of  the 


state  variables  and  their  derivatives  are 
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ay  ■  r  -  r|  -  (Ff  ♦  Lr  ♦  Dr  -  mg)/m,  (2.1) 

%  "  ri  ♦  2rJ  -  (F^  ♦  L+  ♦  Dj/m,  (2.2) 

where  Fy,  Ly,  and  are  the  radial  components  of  the  thrust,  lift,  and 
drag  per  unit  initial  mass  and  g  is  the  gravitational  acceleration  and 
L$»  and  are  the  transverse  components  of  the  thrust,  lift,  and 
drag  per  unit  initial  mass. 

The  transverse  force  is  thus 

(F  -  D)  cos  8  -  L  sin  8, 

and  the  radial  force  is 

(F  -  D)  sin  8  +  L  cos  8  -  mg, 

2  2 

where  g  -  gQrQ/r  ,  gQ  being  the  gravitational  acceleration  at  the  surface 
of  the  planet,  and  rQ  being  the  radius  of  the  planet. 

Let  (r,r,$)  ■  (x^ ,x2 ,x^) ,  the  state  variable  vector,  and  let 

(Pr»B)  ■  (u1,u2),  the  control  variable  vector.  The  state  differential 
equations  are 


dr/dm  ■  r/(mx  ♦  ij)  -  f^x.u,*),  (2.3) 

dr/dm  -  (-gQr2/r2  ♦  r|2  ♦  ((F  -  D)  sin  8  ♦  L  cos  8)/m)/(il  ♦  i2)  . 


f2(x,u,m). 


(2.4) 


d+/dm  »  (-2r#  ♦  ((F  -  D)  cos  $  -  L  sindj/n)/^  ♦  ®2)r  »  f5(x,u,m)j (2.5) 

where  -n^  and  -m2  are  the  mode-1  and  node -2  propellant  flow  rates  per 
unit  initial  nass,  gQ  is  the  gravitational  acceleration  at  the  surface 
of  the  planet,  rQ  is  the  radius  of  the  planet,  and  F,  D,  and  L  are 
thrust,  drag,  and  lift  per  unit  initial  mass.  The  right-hand  sides  of 
these  differential  equations  are  functions  of  the  state  variables,  one 
of  the  control  variables,  the  independent  variable,  and  m^,  m2,  F,  D, 
and  L.  One  must  express  a^,  b^,  F,  D,  and  L  as  functions  of  the  state 
and  control  variables. 

The  thrust  is  given  by 

(y  1)/Y 

F  -  -(2YeRTr/(Ye  -  l)We)1/2(l  -  (pa/pp  *  ®)1/2(i1  ♦  (2.6) 

where  y_  and  W  are  the  specific -heat  ratio  and  the  molecular  weight 

of  the  exhaust,  py  is  the  reactor-exit  stagnation  pressure,  Tr  is  the 

reactor-exit  stagnation  temperature  (given),  p  is  the  atmospheric 

© 

pressure  (a  function  of  altitude),  and  R  is  the. universal  gas  constant.3 
The  specific-heat  ratio  and  the  molecular  weight  of  the  exhaust  are  given 
by 

Y#  ■  (Yj^/CTj  “  *)*!  ♦  T2®2^y2  "  “  X)Kx  * 

i2/(Y2  -  D"2),  (2.7) 

We  *  ("l  *  ♦  *2V» 


(2.8) 
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where  Yj  and  y2  are  the  specific-heat  ratios,  and  and  W2  are  the 
molecular  weights,  of  the  two  propellants  (after  dissociation,  if  it 
occurs.) 

The  propellant  flow  rates  satisfy  the  equations 

.  .  (Y6  ♦  l)/2(r  -  1)  1/2 

»!  ♦  “2  "  “(2/(Ye  ♦  D)  ®  (YeVRV/pA*  (2*9) 

P  -  -(m1/p1  ♦  «2/p2)pr,  (2.10) 

where  is  the  nozzle  throat  area  per  unit  initial  mass  (given  ) ,  P  is 
the  effective  pump  power  per  unit  initial  mass  (given) ,  and  p^  and  p2 
are  the  densities  of  the  two  propellants.4  Substituting  Eq.  (2.7)  and 
(2.8)  into  Eq.  (2.9),  and  eliminating  one  of  the  flow  rates  between  the 
resulting  equation  and  Eq.  (2.10),  one  obtains  an  equation  that  can  be 
solved  numerically  for  the  other  flow  rate,  for  given  reactor  pressure. 
Thus  F  ■  F(pr,pa(r});  py  is  a  control  variable  and  r  is  a  state  variable. 

The  reactor  pressure,  py,  must  satisfy  two  inequality  constraints. 
Since  m^  and  m2  are  both  negative,  py  must  be  bounded  by  the  value  p^ 
for  which  m2  ■  0  (pure  ammonia),  and  the  value  p^  for  which  m^  »  0  (pure 
hydrogen).  To  find  p^,  set  m2  »  0  in  equations  (2.9)  and  (2.10).  After 
eliminating  m^  between  these  equations,  the  resulting  equation  will  be 
solved  for  py.  This  value  will  be  p^.  Similarly,  to  find  pH,  set  ■  0 
in  equations  (2.9)  and  (2.10).  After  eliminating  *2  between  these 
equations,  the  resulting  equation  will  again  be  solved  for  py.  This  time 
this  value  will  be  p^.  Thus  the  constraints  on  the  control  vector  are 


PH  5  Vpr  -  PA* 


(2.11) 
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and 


0  <  Uj-0  <  */2.  (2.12) 

The  drag  and  lift  are  given  by 


D  “  VaM2<:D(M'a)V2'  (2.13) 

L  “  YaPaM2CL(Mfa)Ab/2t  (2.14) 

where  y  is  the  specific-heat  ratio  of  air,  M  is  the  flight  Mach  nunber, 
CD  and  CL  are  the  drag  and  lift  coefficients,  a  is  the  angle  of  attack, 
and  is  the  base  area  per  unit  initial  mass  (given) .  The  Mach  number 
and  the  angle  of  attack  are  given  by 


M  -  (Ha/YaRTa)1/2(r2  ♦  (♦  -  »)2r2)1/2,  (2.15) 

a  *  8  -  tan_1(r/(4  -  «)r),  (2.16) 


where  W  is  the  molecular  weight  of  air,  T  is  the  atmospheric  temper- 
ature  (a  function  of  altitude),  and  u  is  the  planet's  rotational  angular 
velocity.  The  temperature  in  the  i-th  layer  of  the  atmosphere  is  given 
by 

Ta  *  Ti  ♦  (r  -  *1)^,  (2-17) 


and  the  pressure  by 


Pa  *  PitVV 


(*t  *  0), 


(2.1«) 


Pa  -  PA«xp(-(r  -  ri)80wa/RTi)  •  0), 


(2.19) 


where  is  a  constant,  and  and  are  the  temperature  and  pressure 
at  the  base  of  the  layer,  at  r  ■  r^.  Thus  drag  and  lift  are  functions 
of  the  state  variables  and  0,  one  of  the  control  variables. 

Thus  the  right-hand  sides  of  the  state  differential  equations  are 
functions  of  the  state  and  control  variables  and  the  independent  variable 
The  initial  and  final  state  variables  are  given,  i.e. 


*i(*f)  ■  xfi>  -  x0^. 


(2.20) 


where  xn  and  x.  are  specified.  In  the  case  of  a  circular  orbit, 
x  *i 

■  0,  and  ■  (foro^rf^2  w*iere  *£  and  if  ar®  th®  fin*1!  values  of  r 


and  j.  The  final  conditions  oh  the  state  variables  can  be  written  in  the 
form  Gi(x(m^),np  ■  0,  where  G^(x(m),m)  ■  x^(m)  -  x^  for  the  x^  are 
given. 

Having  formulated  a  mathematical  model  of  the  system  and  having 
established  initial  and  final  conditions  on  the  state  variables,  one  must 
define  a  performance  index  which  is  to  be  maximized  or  minimized.  In 
this  case,  the  performance  index  is  I  ■  m^  and  this  performance  index  is 
to  be  maximized. 

Thus  the  problem  is  in  the  form  of  an  optimal  control  problem,  where 
one  maximizes  the  functional  I  ■  subject  to  the  differential  con¬ 
straints  given  by  equations  (2.3)  -  (2.5),  the  inequality  constraints 
given  by  equations  (2.11)  and  (2.12),  and  the  boundary  conditions  given 
by  equations  (2.20). 
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CHAPTER  HI 

THE  OPTIMAL  CONTROL  PROBLEM 

The  control  optimization  problem  is  that  of  finding  a  control  vector 
that  "optimizes"  the  performance  of  the  system.  The  optimal  control 
problem  formulated  in  Chapter  II  is  a  member  of  a  class  of  optimal  con¬ 
trol  problems  considered  below.  In  the  general  problem  discussed  below, 
the  independent  variable  is  chosen  as  t,  for  convenience,  which  is 
equivalent  to  the  mass  variable,  m,  of  the  problem  discussed  in  Chapter 
II. 

Let  the  state  vector  x  be  an  N-vector  and  the  control  vector  u  be 
an  M-vector;  let  x  satisfy  the  set  of  first-order  ordinary  differential 
equations 


dXj/dt  «  f£(x,u,t)  i»l,...,N, 

and  let  the  performance  index  be 

t- 

I  ■  P(x(t£),tf)  ♦  Q(x,u,t)dt. 

I  is  said  to  be  a  functional  of  u.  Let  tQ 
6i(x(tf),tf)  -  0  i-l,...,L  (L  <  N+l) , 


(3.1) 


(3.2) 

and  x(tQ)  be  given,  and  let 

(3.3) 


where  t£  may  or  may  not  be  given.  The  control  vector  u,  nay  have  to 
satisfy  inequality  constraints  as  well.  The  problem  is  to  compute  u(t), 
and  t£  if  it  is  not  given,  such  that  I  is  extrenized. 


This  is  the  Bolza  problea.  Zf  P(x,t)  ■  0,  one  has  the  Lagrange 
problem;  if  Q(x,u,t)  -  0*  one  has  the  Mayer  problem.  The  Bolza  problea 
can  be  put  into  either  the  Lagrange  or  the  Mayer  form,  if  desired.  The 


problea  considered  in  this  thesis  is  a  Mayer  problem,  since  I  ■  t £. 

The  following  are  necessary  conditions  for  the  ext realization  of  the 
4 

performance  index: 

1.  X  is  introduced  as  a  costate  vector  with  N  components,  which 
must  satisfy  the  set  of  first-order  ordinary  differential 
equations 

N 

dX./dt  ■  -3Q(x,u,t)/Sx.  -  l  X,3f, (x,u,t)/3x.  (3.4) 

1  1  j-1  3  3  1 


i**l , « . .  ,N, 


These  equations  are  referred  to  as  the  costate  differential 
equations. 

2.  The  control  vector  u  oust  satisfy  the  set  of  algebraic  equations 

N 

3Q(x,U,t)/3u.  ♦  l  X,3f . (x,u,t)/3u,  *  0  (S.S) 

1  j-1  J  *  1 

1® 1 , . . . ,M, 

if  there  are  no  constraints  on  the  control  vector.  These  equa¬ 
tions  are  referred  to  as  the  optimality  conditions.  Zf  there 
are  constraints  on  the  control  vector,  the  solution  to  these 
equations  may  violate  these  constraints.  Zn  that  case,  one 
chooses  the  control  vector  that  comes  as  close  as  possible  to 
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N 

ox t real zing  the  function  Q(x ,u,t)  ♦  £  A.f  (x,u,t)  without 

j-1  3  3 

violating  the  constraints. 

3.  The  final  conditions 


*i(tf)  ~  9P(x(t£)  itjJ/jXj, 


L 


0 


(3.6) 


i-1 


N, 


mist  be  satisfied,  where  the  q's  are  additional  variables  in¬ 
troduced  which  are  referred  to  as  Lagrange  multipliers. 

4.  If  t£  is  not  given,  an  additional  final  condition 

N 

Q(x(t£)  ,u(t£)  ,t£)  ^(^fiCxftf)  »u(t£)  ,t£)  ♦ 

L 

»P(x(tf)  ,tf)/at  -i^qi»Gl(x(tf)  ,t£)/at  -  0,  (3.7) 

must  be  satisfied.  This  is  referred  to  as  the  transversality 
condition. 

Equations  (3.1)  and  (3.4)  can  be  solved  for  x(t)  and  A(t),  with 
the  control  vector  determined  by  equation  (3.5),  and  with  the  constants 
of  integration,  the  Lagrange  multipliers,  and,  if  it  is  not  given,  t£, 
determined  by  the  initial  conditions  on  the  state  vector,  equations  (3.3), 
and  (3.6),  and,  if  t£  is  not  given,  equation  (3.7). 

The  unknowns  thus  to  be  determined  are  an  N- vector  function  x,  an  11- 
vector  function  X,  an  M-vector  function  u,  L  q's,  N  constants  of  inte¬ 
gration,  and  perhaps  t£.  The  equations  from  which  these  are  to  be  deter¬ 
mined  consist  of  N  state  differential  equations,  N  costate  differential 


equations,  M  optimality  conditions,  L  terminal  constraints,  N  final 
conditions  on  the  costate  variables,  and  if  is  not  given,  a  transver¬ 
sal  ity  condition.  Thus  there  are  the  same  number  of  equations  as  there 
are  unknowns. 


CHAPTER  IV 


THE  NUMERICAL  SOLUTION 

In  general,  the  control  optimization  problem  leads  to  a  nonlinear 
two-point  boundary -value  problem  that  cannot  be  solved  analytically  to 
obtain  the  optimal  control  function. 

There  are  several  methods  for  the  numerical  solution  of  the  Bolza 
problem.  The  purpose  here  will  be  to  present  three  iterative  numerical 
techniques  that  could  be  used  to  solve  the  trajectory  optimization  prob¬ 
lem.  Two  of  the  three  techniques  are  procedures  for  solving  general 
nonlinear  two-point  boundary -value  problems.  The  third  technique  is 
applicable  only  to  control  optimization  problems.  The  methods  are 
briefly  discussed  below  without  the  inequality  constraints  on  the  control 
vector. 

Assuming  that  t^  is  fixed,  x(t^)  is  free,  x(tQ)  is  fixed,  and  tQ  is 
given,  the  two-point  boundary -value  problem  can  be  summarized  by  the 
following  equations: 


dXj/dt  -  ft(x,u,t)  i-l,...,N,  (4.1) 

N 

dX./dt  •  -3Q(x,u,t)/3x.  -  £  X,3f, (x,u,t)/3x.  i*l,...,N»  (4.2) 

1  1  j-1  3  3  1 

N 

3Q(x,u,t)/3u.  ♦  l  X.3f . (x,u,t)/3u.  •  0  i«l,...,M,  (4.3) 

1  j«l  3  3  1 

x(t0)  ■  *0,  (*•*) 

^(tj)  "  ^(^(tpitp/JXj  *  0  ial,..<  ,N.  (4.3) 


The  boundary  conditions  come  from  the  N  initial  conditions  on  x,  given 
by  (4.4),  and  the  N  final  conditions  on  X,  given  by  (4.S).  The  perfor¬ 
mance  index  is 

t- 

I  »  P(x(t£),tf)  ♦  Q(x,u,t)dt. 

Each  of  these  numerical  techniques  are  based  on  the  following  pro¬ 
cedure: 

An  initial  guess  is  used  to  obtain  the  solution  to  a  problem  in 
which  one  or  more  of  the  necessary  conditions  is  not  satisfied.  This 
solution  is  then  used  to  adjust  the  initial  guess  in  an  attempt  to  make 
the  next  solution  come  closer  to  satisfying  all  of  the  necessary  con¬ 
ditions.  If  these  steps  are  repeated  and  the  iterative  procedure  con¬ 
verges,  all  the  necessary  conditions  may  eventually  be  satisfied. 

e 

In  the  steepest  descent  method,  the  control  function  u(t)  is 
guessed.  To  select  a  u(t)  to  begin  the  procedure  ,  one  utilizes  whatever 
physical  insight  he  has  about  the  problem.  Since  the  initial  conditions 
on  the  state  vector  are  given,  the  state  differential  equations  can  be 
numerically  integrated  forward  to  obtain  x(t).  Then  X(t£)  can  be  found 
by  equation  (4.5).  Now  since  the  final  conditions  on  the  costate  vector 
are  known,  the  costate  differential  equations  (4.2),  can  be  numerically 
integrated  backwards  to  obtain  X(t).  The  next  step  is  to  solve  the 
optimality  equations  (4.5),  for  a  new  u(t),  using  the  calculated  x(t)  and 
X(t) .  The  procedure  is  repeated  until  the  difference  between  successive 
anproximstions  becomes  relatively  small  or  until  it  is  evident  that  con¬ 
vergence  will  not  occur,  in  which  case  another  u(t)  needs  to  be  selected. 

The  method  of  steepest  descent  is  generally  characterized  by  ease 
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of  starting;  the  initial  guess  of  u(t)  is  not  usually  crucial.  Yet  as 
a  minimum  is  approached,  the  method  has  a  tendency  to  converge  slowly. 

In  each  iteration,  numerical  integration  of  2N  first-order  differential 
equations  is  required.  In  addition,  the  optimality  equations  must  be 
solved  for  a  suitable  number  of  values  of  t. 

The  variation  of  extremals  technique*  is  an  application  of  the  New- 
ton-Raphson  iteration.  This  technique,  unlike  the  steepest  descent  meth¬ 
od,  can  be  used  to  solve  the  general  two-point  boundary-value  problem. 

In  the  variation  of  extremals  technique,  every  trajectory  generated  by  the 
algorithm  satisfies  equations  (4.1)  -  (4.4),  and  hence  is  an  extremal. 

Here  an  initial  costate  vector  X(tQ)  is  assumed.  Now  the  state  and  co¬ 
state  differential  equations  can  be  numerically  integrated  forwards  to 
solve  for  x(t)  and  X(t),  with  u(t)  being  determined  from  the  optimality 
condition.  The  calculated  X(t^)  is  a  function  of  the  assumed  initial 
costate  vector;  it  does  not  satisfy  equation  (4.5).  Thus  the  set  of  N 
algebraic  equations  represented  by  equation  (4.5)  is  solved  for 
X^tQ) , . . .  ,XN(t0)  by  the  Newton-Raphson  method.  This  set  of  equations 
are  solved  for  the  corrections  to  the  approximate  solution.  The  compu¬ 
tation  of  the  right-hand  sides  of  the  linearized  algebraic  equations  re¬ 
quires  the  numerical  integration  of  the  state  and  the  costate  differential 
equations.  The  computation  of  each  column  of  the  coefficient  matrix  re¬ 
quires  another  integration  of  the  differential  equations.  Thus  at  each 
iteration  tho  2N  differential  equations  must  be  integrated  N+l  times  and 
a  set  of  N  linear  algebraic  equations  must  be  solved. 

In  the  variational  of  extremals  technique,  divergence  may  result 
from  a  poor  guess,  but  once  convergence  begins  (if  it  does),  it  is  rapid. 

Actually,  it  may  be  better  to  guess  x(tp  instead,  since  one  would 


probably  have  more  knowledge  about  the  final  values  of  the  states  from 
the  physical  nature  of  the  problem.  If  one  guesses  x(t^)  instead  of 
A(tg) ,  then  the  costate  and  state  differential  equations  would  now  be 
integrated  backwards  in  time. 

Hie  last  technique  to  be  considered  is  quasilinearization.7  This 
method  can  also  be  used  to  solve  the  general  two-point  boundary-value 
problem.  In  this  method,  an  x(t)  and  A(t)  are  assumed.  The  state  and 
costate  differential  equations  are  linearized  about  this  nominal  tra¬ 
jectory,  bearing  in  mind  that  u,  as  a  function  of  x,  X,  and  t,  can  be 
computed  by  solving  the  optimality  conditions.  One  can  obtain  a  general 
solution  to  the  linearized  state  and  costate  differential  equations  in  the 
following  way.  First,  a  particular  solution  is  obtained  by  numerically 
integrating  the  linearized  equations  with  the  specified  initial  conditions 
on  the  state  variables  and  with  zero  initial  values  for  the  costate 
variables.  N  linearly  independent  solutions  to  the  corresponding  homo¬ 
geneous  differential  equations  are  obtained  by  solving  the  homogeneous 
equations  N  times,  using  a  different  set  of  initial  conditions  each  time. 

A  convenient  choice  of  initial  conditions  would  make  one  of  the  costate 
variables  one,  and  all  of  the  otheT  variables  zero.  The  general  solution 
can  be  written  as  the  sum  of  a  particular  solution  and  a  linear  combin¬ 
ation  of  N  complementary  solutions— N  rather  than  2N,  because  the  par¬ 
ticular  solution  satisfies  the  initial  conditions  on  the  state  variables 
by  itself  and  the  complementary  solutions  all  satisfy  homogeneous  initial 
conditions  on  the  state  variables.  The  coefficients  can  be  determined 
such  as  to  satisfy  the  final  conditions  on  the  costate  variables  by  sol¬ 
ving  a  set  of  linear  algebraic  equations.  The  solution  does  not  satisfy 
the  exact  state  and  costate  differential  equations.  It  satisfies  a  linear 
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approximation  to  those  equations.  This  procedure  is  repeated  using  the 
new  x(t)  and  X(t)  each  time  until  the  difference  between  successive 
approximations  becomes  sufficiently  small  or  until  it  becomes  evident 
that  convergence  will  not  occur,  in  which  case  another  x(t)  and  X(t)  needs 
to  be  selected. 

In  the  quasilinearization  method,  divergence  may  result  fTom  a  poor 
guess.  The  initial  guess  of  the  state-costate  trajectory  is  made  pri¬ 
marily  on  the  basis  of  whatever  physical  information  is  available  about 
the  particular  problem  being  solved.  Each  iteration  requires  the  solu¬ 
tion  of  a  set  of  2N  nonhomogeneous  first-order  linear  differential  equa¬ 
tions,  the  solution  of  a  set  of  2N  homogeneous  linear  first-order  dif¬ 
ferential  equations,  with  N  different  sets  of  initial  conditions,  and  the 
solution  of  a  set  of  N  linear  algebraic  equations. 

The  three  methods  discussed  thus  far  would  all  require  modifications 
to  handle  the  fixed  final-state  problem  where  the  final  value  of  the 
independent  variable  is  not  given  and  also  to  include  inequality  con- 

g 

straints  whenever  needed. 

A  more  flexible  method  is  that  of  finite  differences  which  could 
accomodate  terminal  constraints  of  any  form,  from  the  no-constraint  prob¬ 
lem  to  a  fully  specified  state  for  free  or  fixed  t^.  Replacement  of  the 
derivatives  in  the  differential  equations  to  be  solved  with  finite- 
difference  approximations  results  in  a  set  of  nonlinear  algebraic  equa¬ 
tions  to  be  solved.  These  equations  will  be  solved  with  the  Newton- 
Raphson  iteration.  The  following  discussion  describes  the  procedure. 

Let 


h  -  (tf  -  t0)/K, 


(4.6) 


] 


I  j 
!  '1 


I  I 


where  K  is  a  specified  positive  integer,  and  let 


g3  -  g(tn  ♦  jh) , 


(4.7) 


where  g  *  x,  u,  or  X. 

Using  the  forward-,  backward-,  and  central -difference  approxi- 
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nations 


df(t)/dt  -  (-3f(t)  +  4f(t  +  h)  -  f(t  ♦  2h))/2h  ♦  O0O , 


df(t)/dt  «  (3f(t)  -  4f(t  -  h)  ♦  f(t  -  2h))/2h  +  0(li), 


df(t)/dt  -  (f(t  ♦  h)  -  f(t  -  h))/2h  ♦  0(n  ) , 


(4.8) 


(4.9) 


(4.10) 


the  differential  equations  (3.1)  can  be  approximated  by  the  set  of  alge¬ 
braic  equations 


(xi+l  "  xi’1)/2h  ‘  fi(xJ»u3>Vjh)  "  0  (4*11) 


1*1,..., N  j*l, . . . ,K-1 , 


and  the  differential  equations  (3.4)  can  be  approximated  by  the  set  of 
algebraic  equations 


(-3X°  ♦  A\[  -  X2)/2h  ♦  3Q(x°,u°,t0)/3xi  ♦  yfj  (x°,u°,t0)/3xi*0  (4.12) 


i*l |  •  e  e  fN | 
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(aJ*1  *  Xi"1)/2h  ♦  »Q(x,.u,,t0*jh)/3x1  4 

J1Xk#^(x,*u,*V3h)/ixi  "  0  j-l.---.K-l,  (4.13) 

(3A*  -  4A*”1  ♦  A*"2)/2h  ♦  3Q(xK,uK,t(J+Kh)/3xi  ♦ 

N 

l  Aj3f.(xK,uK,tn+Kh)/3x.  -  0  i*l,...,N,  (4.14) 

j-1  J  3  °  X 

The  central -difference  formula  has  been  used  from  t  between  tQ  and 
tj,  the  forward -difference  foxmula  for  t*t^  and  the  backward-difference 
formula  for  t^t^.  The  equations  represented  by  (4.11)  are  referred  to 
as  the  state  algebraic  equations  and  equations  (4.12)  -  (4.14)  are  re¬ 
ferred  to  as  the  costate  algebraic  equations. 

In  order  to  conform  to  the  same  notation,  equations  (3.3),  (3.5), 
(3.6),  and  (3.7),  can  be  written  as 


G^(xK,tQ+Kh)  ■  0  i*l,...,L, 

N 

3Q(x^,u^,tQ+jh)/3ui  +J^A£3fk(x;*,u*,t0+jh)/3ui  "  0 
i>l,...,M  j«0,...,K, 

N 

aJ  -  3P(xIC,t^Kh)/3xi  ♦^qj3Gj(xK,t0*ICh)/3x1  -  0 
N 

Q(xK,uK,tfi*Kh)  ♦  l  ijf. (xK,uK,t*Kh)  ♦ 
u  i-1  1  1  w 


u  N  ,  w 

3P(x  ,tn+Kh)/3t  -  [  q  3G.(x\t0*Kh)/3t  -  0. 
0  i-1  1 


(4.15) 

(4.16) 


i-l,...,N,  (4.17) 


(4.18) 
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The  Newton-Raphson  method  may  be  used  to  solve  equations  (4.11)  - 
(4.17),  and  if  t£  is  not  given,  equation  (4.18),  for  the  x's,  u's,  X's, 
q's,  and  if  t £  is  not  given,  h.  If  there  are  inequality  constraints  on 
the  control  variables,  one  can  compute  the  trajectory  in  segments.  On 
some  segments,  the  optimality  conditions  do  not  violate  any  of  the 
constraints,  and  may  therefore  be  ignored.  On  other  segments,  some  of 
the  optimality  conditions  cannot  be  satisfied  without  violating  some 
of  the  constraints.  On  these  segments,  the  optimality  conditions  that 
cannot  be  satisfied  without  violating  the  constraints  are  replaced  by 
those  constraints,  treated  as  inequalities . 

A  particularly  important  problem  is  that  of  extremizing  t£  for 
given  final  state.  This  problem  may  be  treated  as  a  Mayer  problem  with 
P*t£,  Q=*0,  and 


Gj (x(t£)  ,t£)  -  x^(t£)  ~  Xj  *  0  i*l,. . .  ,N, 


where  is  given.  In  this  case,  equations  (4.12)  -  (4.14)  become 
xi 

*  jl 

(-3X®  +  4xf  -  X2)/2h  ♦  ^Xj3fj(x°,u°,t0)/3xl  ■  0  1-1,. ...N,  (4.19) 


(xj+1  -  xj_1)/2h  ♦^xjafk(x3,u3,t0+jh)/axi  »  o 


H 


(4.20) 


i-l,...,N  j-1, . . . ,K-1, 

(3X*  -  4X*"1  ♦  X*'2)/2h  ♦^iX|[3fj(xK,uK,t0+Kh)/3xl  -  0 


i*l gee# |N| 


(4.21) 


2 


i  >Jafk(xJ,uj,Vjh)/3u4  •  0  i-l,...,M  j-0,...,K,  (4.22] 

k»l  K  *  0  1 

N 

l  A*f .  (xK.uK,VKh)  ♦1*0.  (4.23] 

i-1  1  1  0 

If  necessary,  central -difference  approximations  can  be  used  for  the 
derivatives  in  equations  (4.19)  -  (4.22).  Equations  (4.11)  and  (4.19)  ■ 
(4.23)  can  be  solved  by  the  Newton-Raphson  method  for  the  x's,  u's,  A's 
and  h.  Again,  if  there  are  inequality  constraints  on  the  control  vari¬ 
ables,  one  can  divide  the  trajectory  into  segments  on  some  of  which  the 
constraints  may  be  ignored  and  on  some  of  which  some  or  all  of  the  op¬ 
timality  conditions  are  replaced  by  constraints. 
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CHAPTER  V 

APPROXIMATE  COMPUTATION  OF  THE  OPTIMAL  TRAJECTORY 

In  order  to  ensure  convergence,  it  may  be  necessary  to  have  a  rea¬ 
sonable  good  initial  approximation  for  the  iterative  procedure  described 
in  the  previous  chapter.  Approximate  trajectories  can  be  obtained  in  the 
following  way. 

It  is  assumed  that  the  effects  of  atmosphere  are  ignored,  i.e., 
p  «  0,  and  the  specific  heat  ratios  of  the  propellants  are  taken  to  be 
equal,  i.e.,  *  y ^  •  y#.  Substituting  equation  (2.9)  into  equation 

(2.6),  one  obtains  F  ■  CjPr,  where 

Cx  -  Y(2/(y  -  l))1/2(2/(y  ♦  l))(Y+1)/2(Y’l)An.  (5.1) 

Thus  the  state  differential  equations  take  the  simpler  form 

dr/dm  -  r/m,  (5.2) 

dr/dm  ■  (-g0rQ/r2  ♦  r$2  ♦  CjPy  sin  S/m)/m,  (5.3) 

di/dm  ■  (-2r$  ♦  Cj,pr  cos  8/m)/mr ,  (5.4) 

where  -m  is  the  total  propellant  flow  rate  per  unit  initial  mass. 

Next  it  is  desired  to  compute  m  as  a  function  of  pr>  Letting 
♦  m2  ■  ra  and  yfl  ■  y,  and  substituting  equation  (2.8)  for  into  (2.11), 
equation  (2.11)  now  is  of  the  form  f(m,m^)  ■  0.  Letting  m^  •  *  -  m^  in 
equation  (2.12)  and  eliminating  ij  between  equations  (2.11)  and  (2.12), 
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and  equating  the  two  of  then,  one  has  an  equation  in  m  only.  This  equa¬ 
tion  is  a  quadratic.  Thus  two  solutions  exist— one  positive  and  one 
negative.  The  negative  solution  is  used  since  n  is  negative  by  defini¬ 
tion.  Thus 

n  -  C3/2C2pr  -  (C2/4C2p2  ♦  C4pJ/C2)1/2,  (5.5) 

where  C 2,  C y  and  are  defined  by 

C2  «  -  P2Wlt  (5.6) 

C3  "  Pplp2CWl  '  V*  (S,7) 

C4  -  y(2/(y  ♦  1))(Y+1)/(Y-1)(P1  -  P2)WlW2An/RTr»  (5,8) 

where  y  is  the  propellant  specific  heat  ratio. 

The  costate  differential  equations  become 

dXj/dm  -  (-X2(2gQrJ/r3  ♦  $*)  ♦  X3(-2r$  ♦  cos  e/n)/r2)/*,  (5.9) 

dX2/dm  •  (-X2  ♦  2X3;/r)/m,  (5.10) 

dX3/dm  »  (-2X2r|  ♦  2Xjr/r)/m.  (5.11) 

The  optimality  conditions  become 


(Xjr  ♦  X2(-gQrJ/r2  ♦  r*2)  -  2X3^/r)(2C5  ♦  (C2  ♦  4C2C4pJ)1/2)  ♦ 


(ICjCjP^/n) (Xj  sin  B  ♦  Xj  cos  B/r)  »  0, 


(5.12) 


tan  B  -  X^r/Xj  »  0.  (5.13) 

B  can  be  obtained  from  Eq.  (5.13).  Equation  (5.12)  can  be  rewritten  as 
a  fourth -order  polynomial  in  pf,  and  this  equation  can  be  solved  itera¬ 
tively  by  the  Bairstow  method,10  which  is  an  application  of  the  Newton- 
Raphson  method;  one  hopes  that  one  at  most  of  these  solutions  will  be 
between  p^  and  p^.  Thus  the  control  variables  can  be  calculated  for 
given  values  of  the  state  and  costate  variables  and  the  independent  vari¬ 
able. 

For  a  set  of  estimated  initial  values  of  the  costate  variables,  the 
state  and  costate  differential  equations  can  be  integrated  numerically. 

A  good  integration  formula  to  use  to  do  this  is  the  Runge-Kutta  method. 
The  right-hand  sides  of  the  differential  equations  (5.9)  -  (5.11)  are 
linear  homogeneous  functions  in  terms  of  X^,  X^,  and  Xj.  Also  the  opti¬ 
mality  condition  (5.12)  is  a  linear  homogeneous  algebraic  equation  in  X^, 
X2»  and  Xy  Equation  (5.13)  is  also  linear  homogeneous  in  X2  and  Xj, 
after  multiplying  it  by  Xj.  Hence,  if  the  initial  values  of  the  costate 
variables  are  multiplied  by  a  common  factor,  the  trajectory  is  unaffected 
and  the  subsequent  values  of  the  costate  variables  are  all  multiplied  by 
the  same  factor.  Therefore  one  can  assign  an  arbitrary  initial  value  to 
one  of  the  costate  variables,  say  Xj.  In  order  to  satisfy  the  trans- 
versality  condition  for  any  desired  final  mass,  one  needs  only  to  multi¬ 
ply  the  estimated  initial  values  of  the  costate  variables  by  the  proper 
factor.  Thus  any  point  on  the  trajectory  can  be  a  final  state.  A  tra¬ 
jectory  computed  in  this  way  is  a  minimum-propellant  trajectory  for  the 


conputed  final  state;  but  this  final  state  is  not  known  beforehand. 

(All  this  is  still  true  if  one  does  not  ignore  the  atmosphere  and  if  one 
does  not  assume  that  the  two  specific  heat  ratios  are  equal;  the  opti¬ 
mality  conditions  and  the  right-hand  sides  of  the  costate  differential 
equations  axe  still  linear  and  homogeneous  in  the  costate  variables.) 

It  follows  from  the  optimality  condition  on  0,  equation  (S.13),  that 
*2  and  Xj  must  always  have  the  same  sign  since  tan  6  >  0  for  every  point 
on  a  realistic  trajectory.  Therefore,  if  X2  and  X^  are  ever  going  to 
change  signs,  they  must  do  it  at  the  same  time,  i.e.  X2Xj  >  0.  One  way 
to  ensure  that  X2  and  Xj  never  have  opposite  signs  is  to  start  X2  and  Xj 
with  the  same  signs  and  ensure  that  they  both  change  in  the  same  direc¬ 
tion  and  that  this  direction  is  away  from  0;  obviously  this  is  not  a 
necessary  condition.  If  X^  changes  in  the  direction  away  from  0,  then 
Xj  and  dXj  should  have  the  same  sign.  Since  dm  is  negative,  one  has 
XjdXj/dm  <0.  By  equation  (5.11),  this  requires  Xj/X^  K  r/*2** 
inequality  implies  the  initial  value  of  tan  B  <  r/rj.  Evidently  at  the 
beginning  of  an  optimized  trajectory,  one  cannot  have  r  »  0,  since  0 
must  always  be  positive.  Therefore,  the  trajectory  must  follow  a  non- 
optimized  boost  phase  which  may  use  either  chemical  or  nuclear  propulsion. 
(One  suspects,  although  this  has  not  been  verified,  that  in  order  to 
achieve  a  circular  orbit,  a  shuttle  would  have  to  execute  a  nonoptimized 
pushover  maneuver.)  Now  to  ensure  that  X2  and  Xj  change  in  the  same 
direction,  one  must  have  (dX2/dm) (dX^/dm)  >  0.  Using  equations  (5.10) 
and  (5.11)  and  assuming  that  X2/Xj  <  r/(r2$),  then  in  order  to  have 
(dX2/dm) (dXj/dm)  >  0,  one  must  have  Xj/X3  <  2$/r.  One  would  also  like  to 
have  the  initial  values  of  Xj  and  X2  chosen  such  that  pr  would  be  within 
its  limits  p„  and  p..  The  optimal  p  decreases  during  ascent.  Therefore, 
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CHAPTER  VI 

NUMERICAL  RESULTS  AND  CONCLUSIONS 

The  purpose  of  this  work  is  to  determine  the  thrust  vector  program 
for  a  mixed-mode  nuclear  shuttle  that  minimizes  the  propellant  required 
to  reach  a  given  orbit.  The  thrust  is  determined  by  the  control  variables 
6,  which  specifies  the  thrust  direction,  and  py,  which  is  the  reactor 
pressure.  Thrust  is  also  a  function  of  reactor  pressure,  reactor  temper¬ 
ature,  the  nozzle  throat  area,  and  the  turbopump  power. 

Application  of  the  methods  of  optimal  control  theory  to  this  problem 
results  in  a  two-point  boundary-value  problem  and  an  associated  algebraic 
problem  as  has  been  shown  in  the  earlier  chapters.  All  of  the  numerical 
data  presented  in  the  following  tables  were  computed  by  the  methods  de¬ 
scribed  in  Chapter  V. 

As  will  be  soon  demonstrated,  the  performance  advantage  of  mixed¬ 
mode  propulsion  turned  out  to  be  quite  small.  Therefore,  it  was  not 
worthwhile  in  this  particular  problem  to  compute  more  accurate  trajec¬ 
tories  using  the  rather  expensive  procedure  described  in  Chapter  IV. 
However,  the  correctness  of  the  computer  program,  0PTC0N ,  written  to 
implement  the  latter  procedure,  was  verified.  In  most  cases  of  course, 
one  will  not  be  so  fortunate  as  to  obtain  a  simple  approximate  solution 
to  an  optimization  problem,  and  the  more  elaborate  procedure  will  have 
to  be  used  from  the  beginning. 

In  Table  1,  the  state  and  control  variables  for  a  minimum -propellant 
trajectory,  computed  by  the  approximate  method  described  in  Chapter  V, 
are  given  as  functions  of  mass.  The  following  set  of  input  parameters 
aro  used: 


SI 


AR  *  10’3  m2/ton, 

P  ■  100  kw/ton, 

Tr  »  3000  K, 

Y  *  1.4. 

For  these  parameters ,  the  reactor  pressures  for  pure  hydrogen  and 
pure  ammonia  flow  are: 

PH  »  6000  KN/m2  pA  -  13,100  KN/m2. 

The  initial  values  for  the  costate  variables  are: 

X1(l)  *  2.27  x  10"11  m-1, 

X2(l)  «  1.69  x  io"7  s/m, 

Xj(l)  «  1  s. 

For  reasons  given  in  Chapter  V,  any  point  on  this  trajectory  may 

be  considered  a  final  state.  For  example,  the  final  mass  for  a  minimum- 

propellant  trajectory  from  radius  r  »  6380  km,  radial  velocity  r  »  0.5 

km/s,  and  angular  velocity  4  ■  7.25  x  10  ”5  s”1  to  r  ■  6590  km,  r  »  1.22 
•  ~4  1 

km/s,  and  4  ■  4.99  *  10  s  ,  is  *  0.6. 

As  Table  1  indicates,  the  initial  value  of  pr  is  a  little  greater 
than  p^,  but  it  drops  to  pH  immediately  and  remains  there  throughout  the 
trajectory.  The  radial  velocity,  r,  decreases  slightly  just  after  igni¬ 
tion  and  then  increases  until  cutoff.  The  angular  velocity,  4,  increases 
until  just  before  cutoff.  Then  it  decreases  slightly. 

Table  2  gives  the  state  variables  at  a  *  .2  as  functions  of  the 
initial  value  of  the  costate  variable  A^,  denoted  by  A ^ Cl) .  All  the  input 
data  are  the  same  as  for  Table  1,  except  for  the  initial  value  of  A^. 

As  Aj  increases,  the  final  values  of  r  and  r  decrease  slightly  and  the 
final  value  of  4  increases  slightly  as  well.  As  A^(l)  increases,  the 
initial  value  of  p  decreases  slightly. 


Table  1.  State  and  Control  Variables  as  Functions  of  Mass  for  a  Typical 
Trajectory 


m 

r,  km 

r,  km  s’1 

Pr,  kN  a'2 

8 

i.o 

6380 

0.500 

0.72Sxl0“4 

7210 

0.823 

0.9 

6420 

0.464 

1.64 

6000 

0.835 

0.8 

6450 

0.491 

2.73 

6000 

0.873 

0.7 

6500 

0.713 

3.87 

6000 

0.935 

0.6 

6590 

1.22 

4.99 

6000 

1.02 

O.S 

6720 

2.11 

5.87 

6000 

1.13 

0.4 

6960 

3.53 

6.67 

6000 

1.24 

0.3 

7350 

5.73 

6.90 

6000 

1.36 

0.2 

7980 

9.18 

6.50 

6000 

1.48 
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Table  2.  State  Variables  at  a  »  .2  as  Functions  of  the  Initial  Value  of 
the  Costate  Variable 


nf1 

r,  tan 

r,  tan  s"1 

♦  .  s-1 

O.lxlO"11 

8000 

9.21 

6.27X10"4 

0.2 

8000 

9.21 

6.28 

0.3 

7990 

9.21 

6.29 

0.4 

7990 

9.21 

6.30 

O.S 

7990 

9.21 

6.31 

0.6 

7990 

9.21 

6.32 

0.7 

7990 

9.20 

6.33 

0.8 

7990 

9.20 

6.34 

0.9 

7990 

9.20 

6.35 

1.0 

7990 

9.20 

6.37 

2.0 

7980 

9.18 

6.47 

t 
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Table  3  gives  the  state  variables  at  a  «  .2  as  a  function  of  the 

initial  value  of  the  costate  variable  X^,  denoted  by  X^Cl).  All  input 

parameters  have  the  same  values  as  those  of  Table  1,  except  for  X^Cl). 

As  1^(1)  increases,  the  final  values  of  r  and  r  increase,  whereas  the 

final  value  of  |  decreases.  The  initial  value  of  py  increases  as  X^Cl) 

increases,  and  for  values  of  X2 Cl)  <  1.4  *  10*^,  the  initial  value  of 

p„  *  p„.  Therefore,  the  shuttle  would  be  running  strictly  on  hydrogen. 

Table  4  shows  the  state  variables  at  m  a  .2  as  a  function  of  A^, 

the  nozzle  throat  area.  In  this  case,  all  the  input  parameters  are  the 

same  as  in  Table  1,  except  for  A  .  The  final  value  of  r  and  4  increase 

n 

as  Ar  increases  while  the  final  value  of  r  first  increases  slightly, 

-3  2 

then  decreases.  For  Ar  <  10  m  /ton,  the  initial  value  of  py  *  pH«  The 

pressure  increases  briefly  after  ignition,  then  drops  back  to  pH»  As 

-3  2 

An  increases  above  10  m  /ton,  the  initial  value  of  py  decreases  and 
-3  2 

for  A  >  1.2  x  10  m  /ton,  this  initial  value  «  pu. 
n  n 

Table  5  gives  the  state  variables  at  m  »  .2  as  functions  of  pump 
power  per  unit  initial  mass,  P.  As  P  increases,  the  final  value  of  r 
first  increases  slightly,  then  decreases.  The  final  value  of  r  increases; 
for  small  values  of  P,  it  increases  very  rapidly.  As  P  decreases  below 
100  kw/ton,  the  initial  reactor  pressure,  py,  decreases,  and  for  P  <  100 
kw/ton,  this  initial  value  of  py  ■  pH»  The  pressure  increases  briefly 
just  after  ignition,  then  falls  back  to  pH«  As  P  increases  above  100 
kw/ton,  the  initial  pr  decreases,  until  tor  P  >  120  kw/ton,  this  initial 
pressure  ■  pH< 

As  is  the  case  for  any  new  idea,  the  feasability  of  mixed-mode  pro¬ 
pulsion  must  be  demonstrated,  both  from  a  theoretical  and  a  practical 
Viewpoint.  Even  though  mixed-mode  propulsion  is  theoretically 
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Table  3.  State  Variables  at  n  >  .2  as  Functions  of  the  Initial  Value  of 
the  Costate  Variable 


*2(1),  S  m"1 

r,  km 

r,  km  s'1 

♦  .  s'1 

1.4xio"7 

7820 

8.95 

7.S2X10”4 

1.5 

7880 

9.04 

7.10 

1.6 

7930 

9.12 

6.77 

Table  4.  State 

Area, 

Variables  at  n  *  .2 

An 

as  Functions  of  the  Nozzle  Throat 

v  -2  ^ 

r,  km 

r,  km  s'1. 

4,  s'1 

0.8x10'" 

7930 

8.62 

5.77X10"4 

0.9 

7980 

8.95 

6.22 

1.0 

7980 

9.18 

6.50 

1.1 

7970 

9.37 

6.74 

1.2 


7960 


9.53 


7.04 
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Tabl«  5.  State  Variables  at  n  *  .2  as  Functions  of  the  Puap  Power.  P 
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advantageous,  it  is  not  a  good  idea  if  it  is  not  practical  to  implement 
or  if  it  does  not  have  substantial  advantages  over  other  propulsion 
systems.  In  designing  a  nuclear  rocket,  the  first  inclination  is  to 
choose  hydrogen  gas  as  the  propellant  because  hydrogen  has  a  low  molecular 
weight— thus  providing  high  exhaust  velocity.  Thus  a  comparison  between 
mixed-mode  propulsion  and  straight  hydrogen  propulsion  is  warranted  here. 

Using  the  same  input  parameters  as  in  Table  1,  but  running  only  on 
hydrogen  propellant,  the  cutoff  radius  v,  is  reduced  from  7980  to  7960 
km,  the  radial  velocity  r,  from  9.18  to  9.17  km/s,  and  the  angular  ve- 
locity  t,  from  6.50  *  io  to  6.37  x  10  s  .  Thus  for  a  particular 
set  of  input  parameters,  it  appears  the  use  of  mixed-mode  propulsion 
only  slightly  improves  the  performance  of  a  shuttle  above  pure  hydrogen 
propulsion.  This  slight  improvement  would  hardly  justify  the  additional 
cost  and  complexity  of  mixed-mode  propulsion. 

A  trajectory  was  computed  for  straight  ammonia  propellant  and  it 
was  found  that  the  performance  for  mixed-mode  propulsion  was  far  better 
than  from  ammonia  propulsion. 

It  must  be  pointed  out  here  that  it  may  be  possible  to  find  input 
parameters  for  which  the  advantage  of  mixed-mode  propulsion  would  be 
greater.  For  example,  if  the  engines  were  very  heavy,  mixed-mode  pro¬ 
pulsion  would  be  more  advantageous.  Chemical  shuttles  are  marginal 
systems,  whereas  nuclear  shuttles  are  not.  Consequently,  nuclear  shuttles 
are  much  less  sensitive  to  changes  in  mission  and  design  parameters. 

This  explains  why  the  gains  from  mixed-mode  propulsion  are  much  larger 
for  chemical  shuttles  than  for  nuclear  shuttles. 

In  order  to  verify  the  correctness  of  0PTC0N,  the  finite -difference 
program,  a  calculation  was  done  for  the  unconstrained  segment  of  the 


trajectory  using  the  same  input  data  as  in  Table  1.  The  program,  OPTCON, 
is  of  course  a  general  one,  and  can  be  used  for  a  wide  variety  of  prob¬ 
lems . 

The  results  from  OPTCON  agreed  closely  with  those  computed  from  the 

method  of  Chapter  V.  Both  methods  were  used  to  compute  a  trajectory  with 

the  initial  state  r  ■  6381  km,  r  ■  0.5  km/s,  4  »  7.249  *  10"5  s”1,  and 

the  final  state  r  ■  6410  km,  r  »  0.4691  km/s,  4  *  1.511  *  10w*  s”*.  The 

final  mass  as  computed  by  OPTCON  after  five  iterations  was  0.9126,  and 

the  final  mass  as  computed  by  the  method  of  Chapter  V  was  0.9125.  Midway 

through  the  OPTCON  trajectory,  pr  was  6,594  KN/m2  and  8  was  0.8182.  Mid- 

2 

way  through  the  other  trajectory,  pr  was  6,565  KN/m  and  8  was  0.8255. 

It  appeared  that  after  five  iterations,  the  OPTCON  solution  was  oscil¬ 
lating  about  a  fixed  point,  but  it  was  not  apparent  that  the  amplitude  of 
the  oscillations  was  decreasing. 


CHAPTER  VII 


RECOMMENDATIONS 

As  was  indicated  previously,  certain  restrictions  were  imposed  upon 
the  problem  from  the  outset,  as  was  appropriate  for  an  initial  trajec¬ 
tory  study.  The  following  extensions  to  the  study  could  be  made,  thus 
making  the  results  more  realistic.  These  studies  will  require  more  effort 
than  the  restricted  problem.  They  are  listed  in  order  of  increasing 
complexity. 

1.  Allowance  for  non -equatorial  launch  sites  and  nonequatorial 
orbits . 

TWo  more  state  variables,  with  the  corresponding  state  and 
costate  differential  equations,  and  two  more  control  variables, 
with  the  corresponding  optimality  conditions,  must  be  introduced. 
The  added  state  variables  are  latitude  and  the  time  derivative 
of  latitude.  Two  additional  angles  must  be  given  in  order  to 
specify  the  thrust  direction. 

2.  Allowance  for  a  coasting  period;  this  is  required  in  order  to 
reach  high  orbits. 

Up  until  now,  it  has  been  assumed  that  the  shuttle  has  been 
operating  under  power  all  the  way.  But  high  orbits  cannot  be 
reached  under  power  all  the  way  because  the  length  of  the  powered 
trajectory  cannot  be  long  enough.  In  order  to  reach  high  orbit, 
one  must  incorporate  a  coasting  phase  in  the  trajectory.  The 
following  describes  how  to  optimize  a  trajectory  that  includes 
a  coasting  phase.  The  first  step  is  to  compute  an  optimal  tra¬ 
jectory  with  some  specified  final  state.  Next,  the  state  vector 
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is  computed  after  a  specified  coasting  time.  Finally,  one  com¬ 
putes  an  optimal  trajectory  from  that  state  to  the  desired  final 
state.  The  final  mass  is  a  function  of  the  state  at  the  begin¬ 
ning  of  the  coasting  period  and  the  length  of  the  coasting 
period.  This  function  can  be  maximized  by  an  application  of  the 
Newton-Raphson  method.  Thus  one  is  faced  with  an  optimal  con¬ 
trol  problem  inside  a  function  optimization  problem. 

3.  Allowance  for  non-ideal -gas  effects  in  the  nozzle  flow;  compu¬ 
tations  should  be  made  for  the  limiting  cases  of  equilibrium 
flow  and  frozen  flow.** 

Before  now,  thrust  has  been  taken  to  be  a  function  of 
reactor  pressure,  reactor  temperature,  the  nozzle  throat  area, 
and  the  turbopump  power,  which  it  is.  But  thrust  is  really  a 
more  complicated  function  than  was  assumed  before.  It  was 
assumed  before  that  the  molecular  weight  of  the  exhaust  gas  was 
a  simple  function  of  the  mixture  ratio  and  independent  of  pres¬ 
sure  and  temperature.  It  was  also  assumed  that  the  specific 
heat  of  the  exhaust  gas  is  independent  of  temperature  and  pres¬ 
sure.  Neither  assumptions  are  correct.  The  molecular  weight 
of  the  exhaust  is  a  function  of  pressure  and  temperature  and 
the  specific  heat  is  a  function  of  temperature  and  pressure. 
Therefore,  the  calculation  of  thrust  as  a  function  of  reactor 
pressure  needs  to  be  changed.  There  are  two  limiting  cases 
which  can  both  be  handled  easily,  and  which  should  both  be  con¬ 
sidered.  They  are  listed  below. 

1.  The  case  of  frozen  flow 

The  composition  of  the  exhaust  gases  is  assumed  to  be 


constant  all  along  the  nozzle  at  the  reactor-exit  equilib 
rium  composition. 

The  case  of  equilibrium  flow 

The  gas  is  assumed  to  have  its  equilibrium  composition  at 
every  point. 
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APPENDIX  A 

PROGRAM  OPTCON 

OPTCON  solves  the  fol lotting  problem:  Let  the  state  vector  x  be  an 
N-vector  and  the  control  vector  u  be  an  M-vector;  let  x  satisfy  the  set 
of  first-order  ordinary  differential  equations 

dXj/dt  -  FtCx,u,t)  i-l,...,N  .  (A. 1) 

Let  tQ,  x(tQ),  and  x(t^)  be  given.  The  problem  is  to  compute  u(t)  and 
t£,  such  that  t^  is  extremized.  This  will  be  referred  to  as  the  extremal 
t£  problem. 

The  following  are  the  necessary  conditions  for  the  extremization 
of  tf: 

1.  The  N-vector  X  must  satisfy  the  set  of  first-order  ordinary 
differential  equations 

N 

dX./dt  -  -  l  X.3F,(x,u,t)/3x.  i»l,...,N,  (A.2) 

1  j-1  3  3  1 

X  is  referred  to  as  the  costate  vector,  and  these  equations  are 
referred  to  as  the  costate  differential  equations. 

2.  The  control  vector  u  mist  satisfy  the  set  of  algebraic  equations 

N 

l  X,3F. (x,u,t)/3u,  ■  0  i»l,...,M,  (A. 3) 

j-1  3  3  1 

These  equations  are  referred  to  as  the  optimality  conditions. 

3.  Since  t.  is  not  given,  the  final  condition 
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N 

i^1Xi(tf)Fi(x(tf)  *u(tf)  'V  ♦  1  *  °* 


(A. 4) 


oust  be  satisfied.  This  is  referred  to  as  the  transversal ity 
condition. 

OPTCON  uses  a  finite  difference  method  to  solve  equations  (A.l)  • 
(A.4).  The  left-hand  sides  of  equations  (A.l)  and  (A. 2)  are  replaced  by 
central  difference  approximations.  Thus  equations  (A.l)  and  (A. 2)  are 
approximated  by  a  set  of  nonlinear  algebraic  equations;  N  state  equations 
at  each  of  K-2  points  and  N  costate  equations  at  each  of  K  points.  In 
addition,  one  has  M  optimality  conditions  at  each  of  X  points  and  a  trans- 
versality  condition.  All  of  these  equations  are  solved  by  the  Newton- 
Raphson  method  for  the  N(K-2)  state  variables,  the  NX  costate  variables, 
the  MX  control  variables,  and  tf. 

Initial  estimates  for  the  state,  costate,  and  control  variables 
must  be  supplied.  To  solve  the  set  of  nonlinear  equations  . . . ,x^) 

>0,...,Fn(x^,...,xn)«0  by  the  Nevrton-Raphson  method,  the  derivatives  of 
the  F's  with  respect  to  the  x's  must  be  calculated.  This  is  done  in 
subroutine  JM  by  use  of  central  difference  formulas.  In  this  case, 
subroutine  G  is  used  to  compute  the  F's.  In  this  case,  the  F's  contain 
derivatives.  Subroutine  G  calls  D  to  calculate  these  derivatives  by 
means  of  finite  difference  formulas.  The  right-hand  sides  of  the  state 
differential  equations  are  calculated  by  DE.  Subroutine  DE  would  need 
to  be  rewritten  for  each  new  optimization  problem;  the  only  changes 
that  would  have  to  be  made  in  the  other  subroutines  are  in  DO  loop  indices 
and  the  like.  NR  calls  subroutine  GE  to  solve  the  linearized  algebraic 
equations  for  the  corrections  to  the  unknowns. 

A  listing  of  the  program  follows. 


IMPLICIT  REAL*8  <A-H,0-Z) 

— K*1+L*(M+N)+(L-2)*N 

—DIMENSION  A(K) ,A1 (K) ,A2<K> ,DACK) ,GA(K,K) ,U<L»M> ,X( (L-2)*N) ,XO(N> 
— XF(N) ,XL(L*N) ,Y(K) ,Y1 (K) ,Y2(K> 

DIMENSION  A( 1 15) ,A1 (1 15) ,A2( 1 15) ,DA( 1 15) ,GA( 115,115), 

1U(30) ,X(39) ,XO(3) ,XF(3) ,XL(45) ,Y( 1 15) ,Y1 (1 15) ,Y2(1 15) 

DIMENSION  P(8) ,R(8) ,T(8) ,TP(8) 

COMMON  TO,XO,XF 

COMMON/COM1/AB,C1 ,C2,C3,C4,G0,0MEGA,P,R,R0,RU,T,TP,XK,XKA,XMA 
EXTERNAL  G 

DATA  TO,XO,XF/1 .03,6.38106,5.02,7.2490-5,6.41006, 

14.69102,1 .5110-4/ 

DATA  AB , AN , GO , OMEG A , RO , RHO 1 ,RH02,RU,TR,W,XK, 

1XKA.XM1 ,XM2,XMA/1 .03,1 .D0,9.807D0,7.272D-5,6.371D6, 

26.821D-1 ,7.D-2,8.314D3,3.D3,1 .05,2*1 .4D0, 8.51600, 
32.01600,2.89601/ 

DATA  P/1 .013D2,2.263D1 ,5.475D0,8.68D-1 ,1  .1090-1 ,5.90-2, 
11.8210-2,1 .0380-3/,  R/6.371D6,6.382D6,6.391D6, 

26. 40306, 6. 41806, 6. 42306, 6. 43206,6. 45D6/,  T/ 

32 . 88 1 502 , 2 . 1 66502 , 2 . 1 66502 , 2 . 2865D2 , 2 .706502 , 

42. 7065D2, 2. 526502,1  .8065D2/,  TP/-6.5D-3,0.D0, 

51 .0-3, 2. 80-3, 0.00, -2.0-3, -4. 0-3,0. DO/ 

DATA  A( 1 )/-6.2500/ 

DATA  XL/2. 270-1 1,1 .690-7,1 .0000,2. 163D-1 1 ,1 .690-7,1 .0000, 
12.0510-11 ,1 .690-7,9.9990-1,1.9350-11,1 .69D-7, 9. 9980-1,1 .815D-11 , 
21 .6910-7,9.9970-1 ,1 .6910-1 1,1 .6910-7, 9. 995D-1 ,1 .5620-1 1 ,1 .6920-7 
39.992D-1 ,1 .429D-I 1 ,1 .6930-7,9.9890-1 ,1 .292D-1 1 ,1 .694D-7, 9. 9850-1 
41  .149D-1 1 ,1 .6950-7,9.9800-1 ,1 .0020-1 1,1  .696D-7,9.975D-1 , 
58.5020-12,1.6970-7,9.9690-1 ,6.931D-12,1 .699D-7, 9. 9620-1 , 
65.3070-12,1  . 7010-7, 9.9530-1.,3.629D-12,1  .703D-7, 9. 9440-1/ 

DATA  U/7 . 20803 ,8.2310-1,7.10303,8. 2330- 1 , 7 .00403 , 8 .2350- 1 , 
16.90903,8.2370-1 ,6 .81 8D3, 8 .2410-1 ,6.73003,8.2450-1 ,6.64603, 
28.2490-1 ,6. 56503,8. 255D-1 ,6.486D3,8.2610-1 ,6.4103,8.2690-1 , 

36 . 33503 , 8 . 27  70- 1 , 6 . 264D3 , 8 . 2860- 1 , 6 . 1 94D3 , 8 . 296D- 1 , 6 . 1 2503 , 

48 . 3080- 1,6.05903, 8. 320-1/ 

DATA  X/6 . 38306 , 4 . 992D2 , 7 . 734D-5 , 6 . 38506 , 4 . 982D2 , 8 . 23 1 D-5 , 

1 6 . 38606 , 4 . 96902 , 8 . 7 380-5 , 6 . 38806 , 4 .9 5402 , 9 . 2580-5 , 6 . 3906 , 4 .93702 
29.7880-5,6.39206,4.91702,1 .0330-4,6.39406,4.89502, 1 .089D-4, 

36 . 39706 ,4. 871 02, 1.1450-4, 6. 39906 ,4.84502,1. 2030-4 , 6 . 40 1 D6 , 
44.81802,1 .2620-4,6.40306,4.78802,1 ,3220-4,6.406D6,4.757D2, 

51 . 384D-4 , 6 . 40806 ,4.72502,1 .4470-4/ 

DO  1  1*1 ,45 

XL( I )*XL( I )/9 .3840-7 

00  2  1-1,45 

A( 1+1 )*XL(  I ) 

00  3  1*1,30 
A( l+46)*U( I ) 

DO  4  1*1,39 
A( l+76)*X( I ) 

CALL  ASS  1 GN  <  6 , 'OPTCON .OUT ' ) 

C1*XK*DSQRT(2./(XK-1 . ) )*(2./(XK-M  . ) )**(0.5*(XK+1 . )/<XK— 1 . ) )*AN 
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C2=RH01 *XM2-RH02*XM1 
C3=W*RH0 1 *RH02* ( XM1 -XM2 ) 

C4=(RH01-RH02)*XM1*XM2*AN»*2*XK*(2./(XK+1.))«» 
1 ( (XK+1 . )/(XK-1 . ) )/(RU*TR) 


C - CALL  NR(G,K,A,DA,GA,A1 ,A2,Y,Y1 ,Y2> 

CALL  NR(G, 1 1 5 ,A,DA,GA,A1 ,A2,Y,Y1 ,Y2) 

WRITE  (6,10) 

WRITE  (6,100)  A(  1 ) 

WRITE  (6,20) 

C - WRITE  (6,110)  (A( I ) , l=2,1+L*N) 

WRITE  (6,110)  ( A( I ) , 1=2,46) 

WRITE  (6,30) 

C - WRITE  (6,120)  (A(I),I=2+L*N,1+L*(M+N)) 

WRITE  (6,120)  (A( I ) , 1=47,76) 

WRITE  (6,40) 

C - WRITE  (6,110)  (A( I ) , la2+L*(M+N) ,K) 

WRITE  (6,110)  (A(l), 1=77, 115) 

10  FORMAT  (1X,1HH) 

20  FORMAT  (1X.1HP) 

30  FORMAT  (1X,1HU) 

40  FORMAT  (1X.1HX) 

100  FORMAT  ( 1 X , 1 P01 0 . 3  > 

C - FORMAT  (N(1X,1PD10.3>) 

110  FORMAT  (3(1X,1P010.3) ) 

C - FORMAT  ( M ( 1 X , 1  PD 10.3) ) 

120  FORMAT  (2(1X,1P010.3>) 

STOP 

END 

SUBROUTINE  NR(F,N,X,DX,FX,X1 ,X2,Y,Y1 ,Y2) 

C - SUBROUTINE  NR  USES  THE  NEWTON-RAPHSON  METHOD  TO  SOLVE  A  SET 

C - OF  N  NONLINEAR  ALGEBRAIC  EQUATIONS  F(X)=0,  WHERE  F  IS  AN 

C - N-VECTOR  FUNCTION  OF  GIVEN  FORM,  AND  X  IS  AN  UNKNOWN  N-VECTOR. 

C - THE  NEWTON-RAPHSON  METHOD  REPLACES  A  SET  OF  NONLINEAR  ALGEBRAIC 

C - EQUATIONS  WITH  LINEAR  EQUATIONS,  THE  SOLUTION  OF  WHICH  SHOULD 

C - BE  CLOSE  TO  THE  ORIGINAL  EQUATIONS.  SUBROUTINE  NR  CALLS  JM 

C - WHICH  USES  A  THREE-POINT  CENTRAL  DIFFERENCE  FORMULA  TO  COMPUTE 

C - THE  COEFFICIENT  MATRIX( JACOBI  AN  MATRIX).  NR  THEN  CALLS 

C - SUBROUTINE  GE  WHICH  USES  THE  GAUSS  ELIMINATION  METHOD  TO  SOLVE 

C - THE  LINEARIZED  EQUATIONS  FOR  THE  CORRECTION  VECTOR.  INITIALLY, 

C - X  MUST  CONTAIN  AN  ESTIMATE  SOLUTION.  THIS  IS  OVERWRITTEN  BY 

C - THE  SUCCESSIVE  APPROXIMATIONS  TO  THE  SOLUTION. 

C -  THE  USER  SUPPLIES  A  SUBPROGRAM  TO  EVALUATE  THE  EQUATIONS 

C - THEMSELVES. 

C -  THE  INPUT  DATA  ARE: 

C - n  NUMBER  OF  EQUATIONS-NUMBER  OF  UNKNOWNS 

C - X  VECTOR  OF  UNKNOWNS 

C -  THE  WORKING  ARRAYS,  WITH  THEIR  DIMENSIONS,  ARE: 

C - OX  N 

C - FX  N  X  N 

C - XI  N 
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C - X2  N 

C - Y  N 

C - Y1  N 

C - Y2  N 

C - DX,  FX,  X,  XI,  X2 ,  Y,  Y1 ,  AND  Y2  MUST  BE  DIMENSIONED  IN  THE 

C - CALLING  PROGRAM. 

C -  F  IS  A  DUMMY  NAME  FOR  A  USER  SUPPLIED  SUBROUTINE  THAT 

C - COMPUTES  F(X) .  THIS  SUBROUTINE'S  ARGUMENT  LIST  IS  (X,Y);  X  IS 

C - AN  N  X  1  INPUT  ARRAY,  AND  Y  IS  AN  N  X  1  OUTPUT  ARRAY  CONTAINING 

C - F(X) . 

C -  THIS  SUBROUTINE'S  ACTUAL  NAME  MUST  APPEAR  IN  AN  EXTERNAL 

C - STATEMENT  IN  THE  CALLING  PROGRAM.  XI ,  X2,  Y1 ,  AND  Y2  ARE  N 

C - VECTORS.  X,  DX,  FX,  XI,  X2,  Y1 ,  AND  Y2  MUST  BE  DIMENSIONED  IN 

C - THE  CALLING  PROGRAM. 

IMPLICIT  REAL *8  (A-H.O-Z) 

DIMENSION  DX(N) ,FX(N,N) ,X(N) ,X1 (N) ,X2(N) ,Y(N) ,Y1 (N), 

1Y2(N) 

EXTERNAL  F 

1  CALL  JM(F,N,X,FX,X1 ,X2,Y1 ,Y2) 

CALL  F(X,Y) 

CALL  GE(FX,Y,N,0X) 

DO  2  1=1 ,N 

2  xm=x<u-oxm 

WRITE  (6,100)  X(1) 

WRITE  (6,110)  (X(l), 1=2,46) 

WRITE  (6,120)  (X(l), 1=47,76) 

WRITE  (6,110)  (X(l), 1=77, 115) 

100  FORMAT  (1X,1PD10.3) 

110  FORMAT  (3(1X,1P010.3)) 

120  FORMAT  (2( 1 X. 1 PD1 0 .3) ) 

DO  3  1=1, N 

IF  (DABS(DX( I )/X( I ) ) .GT.1 .D-6)  GO  TO  1 

3  CONTINUE 
RETURN 
END 

SUBROUTINE  JM(F,N,X,FX,X1 ,X2,Y1 ,Y2) 

C - SUBROUTINE  JM  COMPUTES  THE  DERIVATIVES  OF  THE  N-VECTOR  FUNTION 

C - F(X)  WHERE  X  IS  AN  N-VECTOR.  SUBROUTINE  JM  USES  A  THREE-POINT 

C - CENTRAL -DIFFERENCE  FORMULA  TO  COMPUTE  THE  DERIVATIVES  OF  THE 

C - FUNCTIONS. 

IMPLICIT  REAL *8  (A-H.O-Z) 

DIMENSION  FX(N,N) ,X(N) ,XI (N) ,X2(N),Y1 (N) ,Y2(N> 

DO  1  l»1,N 
X1(I)=X(I) 

1  X2( I )*X(  I ) 

DO  3  J=1,N 
X1( J)=0.999*X( J) 

X2(J)»1.001«X(J) 

CALL  F(X1,Y1) 

CALL  F(X2,Y2) 

00  2  1-1, N 


FX( I , J)a(Y2( I )-Y1 ( I ) )/<0.002«X( J) ) 

X1(J)«X(J) 

X2( J)aX( J) 

RETURN 

END 

SUBROUTINE  GE(A,C,N,X) 

•THE  GAUSS  EL  I  MAT I ON  SUBROUTINE  SOLVES  THE  SET  OF  SIMULTANEOUS 
•LINEAR  ALGEBRAIC  EQUATIONS  AXaC  WHERE  A  IS  A  GIVEN  N  X  N 
•MATRIX,  C  IS  A  GIVEN  N-VECTOR,  AND  X  IS  AN  UNKNOWN  N-VECTOR. 
■A,  C,  AND  X  MUST  BE  DIMENSIONED  IN  THE  CALLING  PROGRAM.  THE 
•USER  MUST  I  INPUT  THE  NUMBER  OF  EQUATIONS  TO  BE  SOLVED.  THIS 
•SUBROUTINE  RETURNS  THE  SOLUTION  OF  THESE  LINEAR  ALGEBRAIC 
•EQUATIONS  TO  THE  MAIN  PROGRAM.  THE  COEFFICIENT  MATRIX  A  AND 
•THE  C  VECTOR  ARE  OVERWRITTEN  DURING  EXECUTION. 

IMPLICIT  REAL *8  <A-H,0-Z> 

DIMENSION  A(N,N) ,C(N) ,X(N) 

M=N-1 

DO  3  K-1 ,M 
L=K+1 

DO  2  I  aL,N 

IF  (DABS(A( I ,K) ) ,LE.DABS( A(K,K) ) )  GO  TO  2 
DO  1  J=K,N 
Al  JaA( I , J) 

AKJ=A(K, J) 

A(  I , J)=AKJ 
A(K, J)=AI J 
CI=C(I) 

CK=C(K) 

C( I )-CK 
C(K)aCI 
CONTINUE 
DO  3  l=L,N 
8aA(l,K)/A(K,K) 

C( I )=C( I )-B*C(K) 

DO  3  J=L,N 

A(  I , J)=A( I , J)-8*A(K, J ) 

X(N)*C(N)/A(N,N) 

DO  5  la1  ,M 
K*N-I 
L=K+1 
X(K)aC(K) 

DO  4  J*L,N 

X(K)=X(K)-A(K,J)*X(J) 

X(K)»X(K)/A(K,K) 

RETURN 

END 

SUBROUTINE  OE(XOOT,XM,U,X) 

■SUBROUTINE  OE  EVALUTATES  THE  RIGHT  HAND  SIDES  OF  THE  STATE 
•DIFFERENTIAL  EQUATIONS  FOR  GIVEN  VALUES  OF  X,  U,  AND  T. 
IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  P(8) ,R(8> ,T(8) ,TP(8) ,U(2) ,X(3) ,XDOT(3) 


SI 


C0MM0N/C0M1/AB,C1 ,C2,C3,C4, 60, OMEGA, P,R,RO,RU,T,TP,XK,XKA,XMA 
PR=U ( 1 > 

BETA=U(2) 

Q=C3/ ( 2 . *C2*PR ) -DSQRT ( (C3/(2.*C2*PR) )**2+C4*PR**2/C2) 

1=0 

IF  (X( 1 ) .LT.R(2) )  1=1 
IF  (I.EQ.1)  GO  TO  3 
IF  (X(1).GT.R<3>)  1=8 
IF  (I .EQ.8)  GO  TO  3 
DO  2  1=2,7 

IF  (X(  1 )  .GT.R(  I )  .AND.Xd  )  .LT.R(  1+1 ) )  GO  TO  3 

2  CONTINUE 

3  TA=T( I )+(X( 1 )-R( I ) )*TP( I > 

IF  (TPd).EQ.O.)  PA=P(I)«DEXP((R(I)-X(1))*G0* 

1XMA/(RU*T( I ) > ) 

IF  (TP(I).NE.O.)  PA=P( I )*(T( I )/TA)**(GO*XMA/(RU* 

1TP( I ) ) ) 

PA=0. 

ALPHA=U(2)-DATAN(X(2 )/(X( 1 )*(X(3)-0MEGA) ) ) 
XNM=DSQRT(XMA*(X(2)**2+(X( 1 )*(X(3)-0MEGA) )**2)/ 

1 (XKA*RU*TA) ) 

FX=C1 *PR-XKA*PA*XNM**2*C0 ( ALPHA , XNM) *AB/2 . 

FY=XKA*PA*XNM**2«CL ( ALPHA , XNM) »AB/2 . 

S8ETA=0SIN(BETA) 

C8ETA=DC0S(BETA) 

XD0T(1 )=X(2)/Q 

XDOT(2)=(-GO*(RO/X( 1 ) )**2+X( 1 )»X(3)**2+(FX*SBETA+ 
1FY*C8ETA)/XM)/0 

XDOT ( 3 ) = ( -2 . »X ( 2 ) «X ( 3 ) + ( FX*CBETA-FY«SBETA ) /XM ) / 

1 (Q*X( 1 ) ) 

RETURN 

END 

SUBROUTINE  D(FU,FX,T,U,X> 

C - SUBROUTINE  D  EVALUATES  THE  DERIVATIVES  IN  EQUATIONS  4.19-4.22 

C - BY  MEANS  OF  CENTRAL  DIFFERENCE  FORMULAS . 

IMPLICIT  REAL *8  (A-H,0-Z> 

C - DIMENSION  FI (N) ,F2(N) ,FU(N,M) ,FX(N,N) ,U(M) ,U1 (M) ,U2(M) 

C - X(N)  ,X1  (N)  ,X2(N) 

DIMENSION  FI (3) ,F2(3) ,FU(3,2) ,FX(3,3) , U(2) ,U1 (2) ,U2(2) , 

1X(3) ,X1 (3) ,X2(3) 

C - DO  1  1=1  ,M 

00  1  1=1,2 
U1(I)=U(I) 

1  U2( I )*U( I ) 

C - DO  2  1=1  ,N 

DO  2  1=1,3 
XI < I )»X< I ) 

2  X2( I )=X( I ) 

C - DO  4  J=1,M 

DO  4  4-1,2 
U1(J)=0.999»U(J) 
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U2(J)-1 .001»U(J) 

CALL  0E(F1,T,U1,X) 

CALL  DE(F2,T,U2,X) 

C - DO  3  l*1,N 

DO  3  1*1,3 

3  FU( I , J)*(F2( I )-F1 ( I ) )/(0.002*U( J) ) 

U1(J)*U(J) 

4  U2(J)*U(J) 

C - 00  6  J*1,N 

DO  6  J“1 ,3 
X1(J)-0.999*X(J) 

X2( J)*1 .001*X(J) 

CALL  DE(F1,T,U,X1) 

CALL  DE(F2,T,U,X2) 

C - DO  5  1*1  ,N 

DO  5  1-1,3 

5  FX( I , J)*(F2( I )-F1 ( I ) )/(0.002*X( J) ) 

X1(J)*X(J) 

6  X2( J)-X(J) 

RETURN 

END 

SUBROUTINE  G(A,Y) 

C - SUBROUTINE  G  COMPUTES  THE  LEFT-HAND  SIDES  OF  THE  STATE, COSTATE, 

C - OPTIMALITY,  AND  TRANSVERSAL  I TY  EQUATIONS.  A  USER-SUPPL I  ED 

C - SUBROUTINE  DE,  WHICH  EVALUATES  THE  RIGHT-HAND  SIDES  OF  THE  STATE 

C - DIFFERENTIAL  EQUATIONS,  IS  REQUIRED.  SUBROUTINE  G  ALSO  USES  ANOTHER 

0 - SUBROUTINE,  D,  TO  CALCULATE  THE  DERIVATIVES  IN  THESE  EQUATIONS. 

IMPLICIT  REAL *8  (A-H,0-Z> 

COMMON  T0,X0,XF 
C - K=1+L*(M+N)+(L-2>»N 

C - DIMENSION  A(K) ,F(N) ,FU(N,M> ,FX(N,N) ,GP(N,L) ,GU(M,L) ,GX(N,L-2) ,P(N,L) , 

C - 1PJ(N),U(M,L),UJ(M),X(N,L),X0(N),XF(N),XJ(N),Y(K) 

DIMENSION  A(1 1 5) ,F(3) ,FU(3,2 ) ,FX(3,3) , 

1GP(3,15) ,GU(2,15> ,GX(3,1 3) ,P(3,15) ,PJ(3) ,U(2,1 5) ,UJ(2) ,X(3, 15) , 

2X0(3) ,XF(3) ,XJ(3) ,Y( 1 1 5) 

H*A(1 ) 

C - DO  1  J*1  ,L 


DO 

1  J 

1,15 

C - DO 

1  I 

1  »N 

DO 

1  1 

1,3 

C— 

— P(  1 

,J> 

A(I+N»(J- 

1)+1) 

1 

P(l 

,J> 

A(I+3*(J- 

1  )+1 ) 

C— 

— DO 

2  J 

1,L 

00 

2  J 

1,15 

C — 

—DO 

2  1 

1  ,M 

00 

2  1 

1,2 

C— 

— U(  1 

,J) 

A(I+M»(J- 

1  )+L»N-M  ) 

2 

U(l 

J) 

AC  M-2*< J- 

1)+46) 

c— 

— DO 

3  1 

1,N 

DO 

3  1 

1,3 

X(l 

.1) 

XO(l) 

C - X(  I  ,L)=XF(  I ) 

3  X( I , 15)*XF( I ) 

C - DO  4  J=2,L-1 

DO  4  J=2,14 

C - DO  4  1*1,  N 

DO  4  1*1,3 

C - x( I , J)=A( l+N*( J-2)+L*<M+N>+1 ) 

4  X(l,J)*A(l+3»(J-2)+76) 

c - oo  13  J*1,L 

DO  13  J*1 ,15 

C - DO  5  1=1,  M 

DO  5  1*1,2 

5  UJ(I)=U(I,J) 

c - 00  6  1*1,  N 

DO  6  1*1,3 
PJ(I)=P(I,J) 

6  XJ(I)*X(I,J) 

TJ*TO+( J-1 )*H 
CALL  DE(F,TJ,UJ,XJ) 

CALL  D(FU,FX,TJ,UJ,XJ) 

c - DO  7  1*1  ,M 

DO  7  1*1,2 
GU(l,J)*0. 

c - oo  7  K=1,N 

DO  7  K=1 ,3 

7  GU(I,J)=GU(I,J)+P(K,J)«FU(K,I) 

IF  (J.EQ.1)  GO  TO  9 

C - |F  (J.EQ.D  GO  TO  11 

IF  (J.EQ.15)  GO  TO  11 
C - 00  8  1*1  ,N 

DO  8  1*1,3  ,  ' 

GX( I , J-1 )*(X( I ,J+1 >-X< 1 ,J-1 ) )/<2.*H)-FC I ) 
GP(I,J)*(PU,J+1)-P(I,J-1))/(2.*H> 

C - 00  8  K*1  ,N 

DO  8  K*1 ,3 

8  GP(I,J>*GP(I,J>+P<K,J)»FX<K,I) 

GO  TO  13 

C - DO  10  1*1  ,N 

9  DO  10  1*1,3 

GP< I ,1 )*(-3.#P< I ,1 )+4.#P( I ,2)-P( I ,3) )/(2.*H) 

C - DO  10  K=1,N 

DO  10  K*1 ,3 

10  GP< I ,1 )*GP( I ,1 )+P(K,1 )#FX(K, I ) 

GO  TO  13 

11  GH-1. 

C - DO  12  1-1  ,N 

DO  12  1*1,3 

GH“GH+P( I , J)*F( I ) 

GP( I , J)*(3.»P( I ,J)-4.*P< I ,J-1 )+P( I , J-2))/(2.*H) 
•00  12  K*1,N 

DO  12  K-1,3 


C 


12  GP(I,J>*GP(I,J>+P<K,J)*FX<K,I> 

13  CONTINUE 
Y( 1 )*GH 

C - 00  15  J-1.L 

00  15  J*1 #15 

c - DO  14  1=1 ,  N 

DO  14  1*1,3 

C - Y(I+N*(J-1)+1>*GP<I,J> 

14  Y(I+3»(J-1)+1)*GP(I,J> 

C - 00  15  1*1  ,M 

DO  15  1*1,2 

C - Y(I+M*(J-1)+L»N+1 )*GU( l,J) 

15  Y(I+2»(J-1)+46)*GU(I,J) 

C - 00  16  J=1  ,L-2 

DO  16  J-1,13 

C - DO  16  1*1, N 

00  16  1*1,3 

C - Y( l+N*( J-1 )+L*(M+N)+1 )*GX( I ,J) 

16  Y( l+3*( J-1 )+76)*GX( I , J) 

RETURN 

END 

FUNCTION  CO ( ALPHA ,XNM) 

IMPLICIT  REAL *8  (A-H,0-Z> 

C0*0. 

RETURN 

ENO 

FUNCTION  CL ( ALPHA, XNM) 

IMPLICIT  REAL*8  (A-H,0-Z> 

a=o. 

RETURN 

END 
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APPENDIX  B 

NEWTON-RAPHSON  METHOD 

Consider  the  set  of  nonlinear  algebraic  equations 


Fj(x1,x2,...,xn)  -  0  j*l»... ,n. 


Let  Cxi»x2»* • • »x*)  be  an  approximate  solution.  The  Taylor  series  ex¬ 
pansion  of  Fj(x2,x2»...,xn)  to  first  order  is 


Fj ^xi »x2 1 


j?  “  Fj^xi»x2»,,,*xn)  *  1  ^xk“xk^aFj^xi*x2,*”»xn^*xk 


j*l>2, . . • ,n. 


If  x*  1 ,x2+* , . . . ,xj+1  is  then  defined  such  that 

VXt'X2 . XJ>  *jtxi*‘  -  . *„)/>**  -  0 


j“l»2,. • . ,n. 


then  Xj+1,  x2+1,...,x*+1  may  be  an  improved  approximation.  Thus  a  set 
of  nonlinear  algebraic  equations  has  been  replaced  by  a  set  of  linear 
equations,  the  solution  of  which  should  be  close  to  the  solution  of  the 
original  equations.  This  set  of  linear  equations  can  be  solved  by  the 
Gauss  elimination  method,  and  the  values  obtained  can  be  used  to  obtain 
another  approximation,  and  so  on,  until  the  difference  between  successive 
approximations  becomes  sufficiently  small,  or  until  it  becomes  clear 
that  convergence  will  not  occur.  In  the  latter  case,  one  starts  over 
with  a  new  initial  estimate. 
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This  procedure  is  called  the  Newton-Raphson  method. 

The  n*n  matrix  having  dF^/dx^  in  the  jth  row  and  the  kth  column  is 
called  the  Jacobian.  It  is  often  necessary  to  calculate  the  Jacobian 
numerically  by  using  central -difference  formulas. 

In  the  case  of  a  single  equation, 

xi+1  »  x1  -  FCxVcdFCxS/dx). 


